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1  Introduction 


1.1  Motion  estimation  methods 

In  the  past,  a  large  number  of  methods  has  been  developed  for  the  estimation  of 
motion  in  image  sequences.  An  early  method  is  to  search  for  the  maximum  of  the 
cross-correlation  function  [1,2],  This  method  might  be  less  suitable  if  the  images 
contain  a  relatively  large  amount  of  low-frequencies  [3].  Other  methods  were  based 
on  a  local  minimization  of  the  displaced  frame  difference  [4, 5],  Block-wise 
minimization  of  a  uniform  error  criterion  was  also  suggested  [3].  The  latter 
block-matching  methods  have  the  advantage  that  they  are  very  suitable  for 
implementation  on  a  single  integrated  circuit.  Therefore,  block-matching  methods 
are  often  preferred  for  video  compression  schemes  [6, 7].  Most  implementations  and 
standards  use  only  one  resolution  level,  but  multiresolution  block-matching  schemes 
are  being  used  more  frequently  now  [8]. 

Some  authors  have  suggested  the  use  of  edge  feature  detectors  [9],  as  well  as  other 
feature  detection  methods  [10].  For  low-resolution  and/or  noisy  images,  these 
methods  might  be  less  suitable.  Linearized  forms  of  the  displaced  frame  difference 
methods  are  the  gradient  methods  [11].  A  difficulty  associated  with  the  gradient 
methods  is  the  aperture  problem,  i.e.  ill-posed  problem  of  finding  the  velocity 
vectors  normal  to  the  object  boundaries.  To  overcome  the  latter  problems,  additional 
smoothness  constraints  have  been  used  [12,  13].  The  aperture  problem  is  not  limited 
to  gradient  methods.  It  also  manifests  itself  in,  for  instance,  block-matching 
schemes. 

Several  methods  try  to  imitate  the  human  visual  system  to  a  certain  extent.  Methods 
using  activation  profiles  were  discussed  by  Waxman  [14].  Spatiotemporal  filters 
using  Gaussian  bandpass  filters  were  shown  to  be  quite  robust  [15],  at  the  cost  of 
high  computational  requirements.  The  local  phase  in  these  filters  can  be  used  to 
arrive  at  stable  motion  estimation  schemes  [16-18].  Combinations  of  correlation 
methods  with  relaxation-labeling  methods  were  suggested  in  [19].  The  trend  is  that 
more  and  more  motion  estimation  schemes  are  first  trying  to  detect  the  contours  or 
objects  and  then  try  to  perform  motion  estimation  [20].  However,  more  work  needs 
to  be  done  to  arrive  at  stable  and  efficient  methods. 


1.2  Desired  properties 

Desirable  properties  of  methods  which  are  intended  to  determine  motion  between 
images  are  described  by  [21].  The  methods  which  are  pertinent  to  our  problems  are. 
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•  Multiresolution  property.  The  motion  estimation  methods  should  be  able  to 
estimate  velocity  vectors  over  a  large  magnitude  range.  In  addition,  the 
methods  should  be  able  to  handle  a  large  range  of  object  sizes. 

•  Modest  memory  requirements.  Global  optimization  methods  often  lead  to  the 
solution  of  a  large  number  of  linear  equations.  The  solution  to  this  system  of 
equations  is  only  possible  for  a  small  number  of  unknowns  if  direct  matrix 
inversion  methods  are  used  due  to  the  prohibitive  storage  requirements. 

•  Computational  efficiency.  The  real-time  constraint  for  contemporary  video 
formats  imposes  considerable  difficulties,  even  for  relatively  simple  methods. 

•  Subpixel  accuracy.  Most  hardware  implementations  only  calculate  integer 
displacements,  which  might  not  be  sufficient. 

•  Convergence.  Because  of  the  large  number  of  unknowns,  many  motion 
estimation  schemes  are  based  on  iterative  schemes.  The  convergence  of  these 
schemes  is  not  always  garanteed,  however.  This  might  lead  to  strange  motion 
vector  estimates. 

•  Robustness.  The  main  application  area  of  our  motion  estimation  algorithms  is 
motion  detection.  The  primary  objective  is  to  create  a  robust  system,  rather 
than  to  provide  a  system  which  yields,  for  example,  very  accurate  motion 
vectors. 
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2  Motion  estimation  with  global  minimization  schemes 


2,1  Iterative  methods  based  on  global  error  minimization 

Most  of  the  motion  estimation  schemes  are  based  on  the  formulation  of  some  kind 
of  continuity  constraint,  such  as  the  continuity  of  image  intensity.  A  property  of  the 
continuity  constraints  is  that  they  impose  local  requirements.  A  disadvantage  of 
methods  using  local  requirements  is  that  they  are  extremely  sensitive  to  noise.  In 
this  section  we  try  to  optimize  the  continuity  constraints  over  the  complete  image. 
An  extension  to  incorporate  continuity  with  respect  to  the  time  coordinate  is  also 
rriven  The  methods  amount  to  the  solution  of  a  large  system  of  linear  equations. 
Iterative  techniques  are  essential  for  the  solution  of  the  problem,  because  the 
required  computer  storage  and  computation  time  of  direct  matrix  inversion  methods 
are  prohibitive.  The  solution  of  the  system  of  equations  will  be  obtained  with  the  aid 
of  a  continuous  version  of  the  conjugate  gradient  technique.  The  latter  technique 
converges  to  the  global  optimum  for  the  pertaining  quadratic  error  criterion. 

For  other  error  criteria,  the  conjugate  gradient  technique  might  be  less  suitable. 
However,  the  quadratic  error  criterion  is  a  natural  choice  for  most  physical  systems 
because  it  minimizes  the  energy  difference  between  two  quantities.  If  the  data 
suffers  from  severe  outliers,  then  the  robustness  might  be  increased  with,  for 
example,  the  mean-absolute-value  error  criterion.  In  this  report  only  the  widely  used 
mean-squared  error  criterion  is  used  because  it  leads  to  sufficiently  accurate  and 
stable  results  in  nearly  all  cases  of  practical  interest. 


2.2  Vectorial  displaced  frame  difference 

The  detection  of  motion  can  be  based  on  Displaced  Frame  Difference  (DFD) 
methods.  Picture  element  (Pel)-recursive  DFD  schemes  have  been  derived  in  the 
past.  These  schemes  try  to  estimate  the  image  motion  on  a  pixel-by-pixel  basis  by 
minimizing  a  local  error  criterion,  the  Displaced  Frame  Difference  [4].  A 
disadvantage  is  the  poor  convergence,  if  convergence  is  obtained  at  all.  An 
additional  disadvantage  is  the  susceptibility  to  noise.  In  this  section  we  try  to  amve 
at  a  converging  iterative  scheme  with  improved  stability  where  we  minimize  a 
vectorial  version  of  the  DFD-error  in  a  least-mean  square  sense,  i.e.,  globally,  with  a 
conjugate  gradient  technique.  The  only  requirement  is  that  the  DFD  can  be 
linearized  with  respect  to  the  pixel  displacement. 
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2.2.1  Derivation 

The  starting  point  of  our  derivation  is  the  formulation  of  continuity  of  intensity 
according  to 

/(x,f)  =/(x  +  d,t  + Af).  (2.1) 

The  deviation  from  this  equation  is  expressed  in  the  Displaced  Frame  Difference  [4]: 

dfd(x,  d)  =  /(x,  f)  —  /(x  +  d,  i  +  Af) .  (2.2) 

To  arrive  at  the  desired  results  we  introduce  a  vectorial  version  D  of  the  standard 
dfd  according  to 

Ti  =  {D,,Dy),  (2.3) 

with 

{Dj:,  Dy}  =  dfd(x,  {do:,  dy}),  (2.4) 

in  which 

d,  =  (4,0)  (2.5) 

dy  =  (0,4).  (2.6) 

The  value  of  the  displacement  d  will  be  searched  for  in  an  iterative  way  [22].  The 
approximations  of  d  at  iteration  n  will  be  denoted  by  d(^l .  A  global  error  criterion 
at  iteration  n  is  introduced  as 


where  F(")  =  D(x,  d^”)).  If  =  0,  we  have  =  0,  Vx  G  S,  and  we  have 
solved  our  problem,  also  in  a  local  sense.  We  assume  that  D  can  be  linearized,  i.e., 

D(x,  adi  + /3d2)  =  q;D(x,  di)  + /3D(x,  d2),  Vdi,d2.  (2.8) 

This  assumption  can  be  justified  on  grounds  of  proportionality  of  the  dfd  with 
respect  to  image  shifts  in  a  certain  direction.  The  displacement  field  d("^  is 
expanded  in  vectorial  Fourier  components  =  <j!)("^(k)  according  to 

/OO 

(2.9) 

-CXO 

with  k  =  (4,  ky)  the  spatial  frequencies  in  the  (x,  y)-directions.  The  constraint  that 
d("^  is  real  leads  to  symmetry  properties  in  the  spectral  domain.  In  going  from  the 
(n  —  l)-st  step  to  the  n-th,  we  take 


^(n)  ^  +^(")g(n) 


(2.10) 
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in  which  are  the  search  directions  and  is  a  variational  parameter  to 
minimize  We  obtain 

p(n)  _  p(n-l)  _ 


fW  =  _D(^x,y  . 

The  error  criterion  can  be  written  as 
j(n)  ^  f 

J  x£S 

with 

^(n)  ^ 


BW  =  f  IlfHfdx. 

ax65 


The  minimum  of  as  a  function  of  r]^^\  is  obtained  at 


Using  this  particular  value  of  we  get 

j(ri)  _  j(n-l) 


(^(n))2 

5(")  ’ 


aM 

p(n)  _  p(n-l) _ f(n) 

The  latter  expression  leads  to  the  orthogonality  property 

[  f(")  ■  =  0, 


which  will  be  used  later.  Using  the  expression  for  and  interchanging 


integrations,  we  can  write 


s(n-i)*  .g('^)dk, 
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in  which 


s 


(")  =  _ 


(2.22) 


The  orthogonality  property  for  and  leads  to  an  orthogonality  property  in 
the  spectral  domain: 


•  =  0. 


(2.23) 


Up  til  now,  the  search  directions  have  been  completely  arbitrary.  We  take  the 
conjugate  gradient  directions 

+ (2.24) 


while 

g(i)  =  s(0). 


(2.25) 


Substituting  these  expressions  in  Eq.  (2.21)  and  using  the  orthogonality  relation  Eq. 
(2.20),  we  obtain 


^(n) 


(2.26) 


A  preconditioned  version  of  the  conjugate  gradient  scheme  is  obtained  by 
introducing  a  preconditioner  P  and  its  adjoint  P*  [23]  according  to 

„(n)  ^  p*pg(n-l)  >  2,  (2.27) 


while 


g(l)  ^  p*ps(0). 


(2.28) 


The  preconditioner  is  chosen  in  such  a  way  that  it  approximates  the  inverse  operator 
of  the  problem  at  hand.  Substituting  these  expressions  in  Eq.  (2.21)  and  using  the 
orthogonality  relation  Eq.  (2.20),  we  obtain 

/OO 

||Ps("-i)||-dk.  (2.29) 

-OO 

Useful  results  are  often  obtained  by  using  the  inverse  of  the  least-squares  diagonal 
[24]. 
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2.2.2  Properties 

A  few  observations  can  be  made  regarding  the  properties  and  the  implementation  of 
the  iteration  scheme. 

•  Provided  the  linearity  assumption  holds,  the  error  decreases  at  each 

iteration  step  provided  >  0.  From  Eq.  (2.29)  we  can  see  that  is  real 
and  positive,  unless  vanishes.  However,  in  the  latter  case  we  have 

arrived  at  the  exact  minimum  in  the  iteration  n  —  1 . 

•  The  evaluation  of  D  for  different  displacements  can  be  performed  with 
interpolation  of  7,  such  as  with  bilinear  interpolation. 

•  The  major  computational  task  at  each  iteration  is  concerned  with  the 
evaluation  of  four  two-dimensional  Fourier  integrals  which  can  be  computed 
efficiently  with  forward  and  inverse  FFTs. 

•  Memory  requirements  are  modest:  the  iteration  scheme  is  a  full  2D 
least-squares  minimization  scheme  without  the  memory  requirements  of 
standard  matrix  inversion  methods. 

•  Motion  compensation  is  possible  by  using  the  results  for  the  kx  =  0.,ky  —  0 
component. 

•  The  problem  should  be  scaled  properly  because  of  the  non-linear  dependence 
of  D  on  X. 

2.2.3  Iteration  scheme 

The  complete  iteration  scheme  can  be  found  below.  The  iteration  scheme  will  be 
started  with  zero  estimates. 

•  Initialization 

=  0 

=  D(x,0) 

j(0)  ^  f  ||f(0)||2^x 

J  x€S 

=  0 
n  =  0 


n  =  n  +  1 


•  Start  iteration 
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p(”) 

j(") 

•  End  iteration 

d(") 


2.3  Gradient  method  with 

The  detection  of  motion  can  also  be  based  on  spatial  and  temporal  image  gradients 
[11].  The  equivalence  of  gradient  methods  and  correlation  based  methods  has  been 
noted  [25].  The  gradient  schemes  try  to  estimate  the  pixel  motion  by  minimizing  a 
local  error  criterion  based  on  the  continuity  of  image  intensity  [11].  The  velocity 
vector  is  obtained  by  a  global  minimization  of  the  gradient  equation  with  a  conjugate 
gradient  technique.  Smoothness  constraints  are  imposed  which  simultaneously 
minimize  the  variation  in  the  estimated  velocity  in  space  and  time.  Additionally, 
simplified  schemes  are  given  which  are  based  on  the  normal  component  of  velocity. 


,(")  =  ^ 


=  0 


smoothness  constraints 


2.3.1  Derivation 

We  start  our  derivation  by  recalling  the  so-called  gradient  equation,  which  can  be 
obtained  from  linearization  of  Eq.  (2.2)  (see  [11]): 

at/(x,  t)  +  v(x,  t)  •  VJ(x,  t)  =  0,  (2.30) 

where  we  have  used  v  =  df  and  dt  denotes  differentiation  with  respect  to  time. 
Pixel-wise  solution  of  this  equation  for  v  leads  to  serious  problems:  if  the  gradient 
approaches  zero  then  the  result  becomes  unbounded.  To  overcome  the  latter 
problems,  we  define  a  global  error  criterion  which  is  minimized  in  a  least-squares 
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sense  by  imposing  smoothness  constraints,  as  in  the  error  criterion  of  [12].  In  this 
case,  however,  we  also  impose  a  smoothness  constraint  in  the  time  direction.  The 
total  error  criterion  becomes; 

J=  /  dt  f  +  + 

JteT  ^ 

a[\\Vv,f  +  \\Vvyf-]+P\\dtvf}dx,  (2.31) 


The  error  criterion  is  minimized  with  a  conjugate  gradient  technique.  The  global 
error  criterion  at  iteration  n  is  written  as 

=  [  dt  [  +  ||Fi^)l|2  +  ||f(")112  +  IIFj'^f )  dx, 

JteT  Jxes  ^  (2.32) 

where 

=  dtl{x,  t)  +  (x,  f)  ■  V/(x,  t),  (2.33) 


Fg  = 


(2.34) 


f(^)  =  (2.35) 

and  =  (ui"\  4"^)  the  approximation  of  v  at  iteration  n.  The  velocity  field 
is  expanded  in  vectorial  Fourier  components  0^"^  =  (/)^"^(k,  w)  according  to 


/oo  roo 

du  xeS,t€T,  (2.36) 

'CO  —CO 


with  k  —  {kj;,  ky)  the  spatial  frequencies  in  the  {x,  y)-directions.  The  constraint  that 
is  real  leads  to  symmetry  properties  in  the  spectral  domain.  In  going  from  the 
(n  —  l)-st  step  to  the  n-th,  we  take 

0(n)  ^  ^(n-l)  ^^(n)g{n)_  (2.37) 


We  obtain 


with 


{fW.rW  J  =  (2.38) 


/OO  roo 

•OO  j  —CO 


(2.39) 


f(") 

^,y 


x,y 


)gj(k-x- 


(2.40) 
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f(") 


dujuj 


(2.41) 


The  error  criterion  can  be  written  as 

j(n)  ^  j(n-i)  _  (2.42) 

with 

j^{n)  ^  f  ( p{n-l)j{n)  ^  p(n-l)  .  f(.)  +  p(n-l)  .  f(n)  ^  p(n-l)  ,  f(n)\ 

iter  2x65  ^  (2.43) 


and 


/  / 

h&T  Jxes 


(|/Wi2  +  ||fW||'  +  llfi")|i'  +  ||fi^"^ll") 


dx. 


The  minimum  of  as  a  function  of  is  obtained  at 

A{n) 


(n)  _ 


J5H' 


Using  this  particular  value  of  we  get 


j(n)  _  j(n-t) 


SH 


(2.44) 


(2.45) 


(2.46) 


and 


^(r 


S(") 

The  latter  expression  leads  to  the  orthogonality  property 


[  dt  [  dx  = 

JteT  Jx€S  ^  ^ 


(2.47) 


0, 

(2.48) 


which  will  be  used  later.  Using  the  expression  for  and  interchanging 
integrations,  we  can  write 

/oo  rco 

duj  (2.49) 

-OO  J  —OO 


in  which 


dt  [  -  ia^l'^k  ■  F^  + 

2x65  ' 


(2:50) 
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s(”)  — 


dt 

't£T  Jxes 


(2:51) 


with  Ft  =  {Ftx,  Fty).  The  orthogonality  property  for  and 

{/("),  J  leads  to  an  orthogonality  property  in  the  spectral  domain; 


du 


d.n)*  .  g(n)^j^  ^  0. 


(2.52) 


Up  til  now,  the  search  directions  g("^  have  been  completely  arbitrary.  We  take  the 
conjugate  gradient  directions 


g 


(n)  _  g(n-l)  ^ 


>  2, 


(2.53) 


while 

g(i)  =s(°). 


(2.54) 


Substituting  these  expressions  in  Eq.  (2.49)  and  using  the  orthogonality  relation  Eq. 
(2.48),  we  obtain 

/oo  rOQ 

du  /  l|s("-i)fdk.  (2.55) 

-OO  j  —OO 

A  preconditioned  version  of  the  conjugate  gradient  scheme  is  obtained  by 
introducing  a  preconditioner  P  according  to  [23] 

g(n)  ^  p*ps(rr-l)  p  >  2,  (2.56) 


while 


g(t)  ^ 


(2.57) 


where  P*  is  the  adjoint  of  P.  Substituting  these  expressions  in  Eq.  (2.21)  and  using 
the  orthogonality  relation  Eq.  (2.20),  we  obtain 

/CO  POO 

du  /  !lPs(”-^)fdk.  (2.58) 

-CO  J  —oo 


Useful  results  are  often  obtained  by  simply  using  the  inverse  of  the  least-squares 
matrix  diagonal  [24],  such  as  in 


dt 


UteT 


f  (|9x,y2p -hoUkf -h/3u;^)  dx 

7xe5  J 


-1/2 


.(n) 

x,y’ 


(2.59) 
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2.3.2  Properties 

•  The  error  decreases  at  each  iteration  step  provided  .4^”^  >  0.  From  Eq. 

(2.58)  we  can  see  that  is  real  and  positive,  unless  vanishes. 

However,  in  the  latter  case  we  have  arrived  at  the  exact  minimum  in  the 
iteration  n  —  1. 

•  The  major  computational  task  at  each  iteration  is  concerned  with  the 
evaluation  of  four  two-dimensional  Fourier  integrals  which  can  be  computed 
efficiently  with  forward  and  inverse  FFTs. 

•  Memory  requirements  are  modest;  the  iteration  scheme  is  a  full  2D 
least-squares  minimization  scheme  without  the  memory  requirements  of 
standard  matrix  inversion  methods. 

•  Motion  compensation  is  possible  by  using  the  results  for  the  —  0,ky  =  0 
component. 

•  It  should  be  noticed  that  v  is  a  vector  with  two  spatial  components.  Therefore, 
the  total  number  of  unknowns  is  twice  the  number  of  equations  if  standard 
FFTs  are  used  for  computing  the  spatial  and  spectral  integrals.  As  a 
consequence,  the  number  of  spectral  components  should  be  limited  to  a 
maximum  of  half  the  number  of  v-vectors. 

•  For  subsequent  frames,  the  iterative  scheme  can  be  started  with  the  result  v  of 
the  previous  frame. 

2.3.3  Iteration  scheme 

The  starting  value  of  the  iteration  scheme  is  not  critical  because  the  convergence  is 
assured.  For  simplicity,  the  iteration  scheme  will  be  started  with  zero  estimates. 

•  Initialization 

=  0 
=  dtl 
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Start  iteration 


— 


/ 

'ter  Jxes 

j  dtj 

IteT  Jxes 


-F^n)dj  _  i«t/2k  .  f(")  + 


g 


^(n) 
{n)  _  / 


n  =  n  +  1 

f^OO  /’CO 

du  /  ||Ps("-')||2dk 

-oo  «/  — oo 


p2g(0) 


n  =  1 


F2s(^-i)  +  n  >  1 


/OO  roo 

du)  / 

-OO  J  — oo 

/oo  /•oo 

du;  /  k3We'('^"-^‘’dk, 

-OO  «/  — oo 
/»oo  roo 

dujio  gHe'(>^-^-“')dk, 
7—00  — OO 


p(n)  = 


dt  [  +  llf'"^ll'  +  llfi"^ll-  +  llft^”^ip) 

^(n) 


'ter  dxe5 


dx 


M  _ - 


0(n)  ^  +^WgH 

pin)  _  p(n-l)  _  ^{n)  j{n) 

^  x,y,t  —  ^  x,y,t  ^  ^x,y,t 
j(n)  ^  j(n-t)  _  7y(")yl(") 


End  iteration 


r(") 


r»oo  roo 

du  / 

J  —OO 


-OO  J  — oo 


2.3.4  Results 

In  this  section  we  show  an  example  of  the  algorithm  of  subsection  2.3.3  for  two 
images.  Consequently,  we  set  P  =  0.  The  velocity  estimates  are  obtained  after 
low-pass  filtering  of  the  images.  It  can  be  seen  in  Figs.  2.2  and  2.3  that  the  estimated 
velocity  of  the  middle  part  of  the  car  is  different  from  the  estimated  velocity  near  the 
boundaries.  Only  the  object  boundaries  normal  to  movement  are  detected.  This  is  a 
well-known  problem  [11]  which  can  be  attacked  by  incorporating  other  methods, 
such  as  contour  detection.  For  some  applications,  such  as  motion  detection,  it  is  not 
very  important  to  know  the  motion  vectors  that  are  parallel  to  the  object  boundaries. 
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Figure  2.1:  Two  consecutive  image  frames. 


Figure  2.2:  Computed  velocity  of  Fig.  2.1 
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Figure  2.4:  Image  of  Fig.  2.1  with  superimposed  and  thresholded  velocity. 


2.4  Scalar  velocity  method 

It  is  possible  to  reduce  the  number  of  unknowns  by  reducing  the  vectorial  velocity  to 
the  scalar  velocity  component  in  the  direction  of  the  spatial  gradient  (see  also:  [11]). 
Without  using  the  smoothness  constraints,  the  continuity  requirement  of  the 
previous  section  can  be  written  as 

atl(x,i)  +  ||v(x)l|  ||V/(x,f)||cos0(x)  =  0,  (2.60) 

with  0  =  Z{\,  V/)  the  angle  between  the  gradient  and  the  velocity  vector.  If  we 
write 


^(x)  =  l|v(x)||  cos  i9(x),  (2.61) 

we  obtain 

dtl{x,t)  +t;(x)||V/(x,f)l|  =  0.  (2.62) 

Instead  of  solving  for  the  two  vectorial  components  of  v  we  now  try  to  solve  for  the 
scalar  component  v.  The  iteration  scheme  will  be  given  without  derivation. 

2.4.1  Iteration  scheme 


•  Initialization 


=  0 
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=  dtl 

j(0)  ^  f  |ir(0)|2^x 
Jxes 

=  0 

n  =  0 


•  Start  iteration 


sW  =  -/ 

Jx€S 

n  —  n+  1 

/OO 

-OO 


fin) 


A(^) 

s(") 


=0(^-1) +77 

pin)  _  ^(n-l)  _  ^(n) y(n) 
j(n)  _  j(n-l)  _  ^{n)j^in) 


•  End  iteration 


= 


e^^-^dk 


2.4.2  Example 

The  gradient  method  based  on  the  scalar  velocity  was  applied  to  recorded  image 
frames.  Two  consecutive  frames  can  be  found  in  Fig.  2.5.  In  Fig.  2.6,  the  computed 
velocity  is  shown.  Figure  2.7  shows  the  first  image  frame  in  which  a  detection  limit 
has  been  used  to  indicate  image  motion.  The  images  have  a  size  of  256  x  256 
pixels.  About  six  iterations  were  used.  This  is  not  critical,  however:  often  only  one 
iteration  leads  to  sufficient  convergence. 


TNO  report 


FEL-96-B265 


21 


Figure  2.6:  Normal  velocity  computed  from  the  two  consecutive  image  frames. 
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Figure  2.7:  First  image  frame  with  indicated  image  velocity. 

2.5  Scalar  velocity  method  with  smoothness  constraint  in  time 

As  a  last  example,  a  scheme  is  given  which  uses  scalar  velocity,  Fourier  descriptions 
in  time  only,  and  a  smoothness  constraint  in  time.  Then,  the  error  criterion  becomes 

J=  /  dt  [  \[dtI  +  v\\yi\\]'^  +  0\dtv\‘^}dK.  (2.63) 

Jt&T  Vx€5  ^ 

The  iteration  scheme  will  be  given  without  derivation.  An  advantage  of  this  scheme 
is  the  relatively  low  computational  complexity,  because  the  temporal  Fourier 
transforms  are  one-dimensional  and  relatively  short.  For  very  short  FFT’s,  the 
pertinent  matrix  may  consist  of  only  ones  and  zeros,  which  makes  the  scheme 
interesting  for  hardware  implementations. 


2.5.1  Iteration  scheme 


•  Initialization 


=  0 
=  dtl 
=  0 


dx 


,^(0)  ==  0 
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n  =  0 


•  Start  iteration 


^(n)  ^  ^(n-1)  +^(n)^(n) 
p{n)  _  pin-l)  _  ^{n)j{n) 
pin)  _  ^^(’^-1)  _  ^(n)  ^in) 
j(n)  _  j(n-l)  _  r]^n)^{n) 


•  End  iteration 


e-iujt 


duj 


2.5.2  Multi-resolution  analysis 

The  continuity  equations  of  the  preceding  section  are  based  on  linearization  of  the 
displaced  frame  difference.  Therefore,  large-scale  movements  can  only  be 
determined  from  low-pass  filtered  images.  In  the  present  section,  we  define  a 
multi-resolution  scheme  where  image  motion  is  determined  from  subsequent 
application  of  one  of  the  gradient-based  conjugate  gradient  schemes  to  low-pass 
filtered  images.  Each  subsequent  conjugate  gradient  calculation  is  performed  with 
increasing  cut-off  frequency  fc,  i.e.,  on  a  finer  scale.  After  application  of  the 
conjugate  gradient  scheme,  we  fall  back  on  the  displaced  frame  difference  by 
performing  interpolation  in  the  original  image.  The  interpolated  image,  which  is 
based  on  the  velocity  vector  defined  in  the  previous  step,  is  used  as  a  basis  for 
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determining  image  motion  on  a  finer  scale.  The  basic  idea  is  that  we  try  to  estimate 
the  motion  vector  needed  to  get  from  the  first  image  to  the  second  image.  The 
estimated  second  image  is  denoted  by  using  a  tilde  above  the  different  symbols.  The 
procedure  can  be  summarized  in  the  following  scheme  (the  symbol  **  denotes  a 
two-dimensional  convolution): 

•  initialization 

n  =  0 
==  0 
+  At)  = 

=  fcMn 

•  Start  iteration 

t  +  At)  ^  LPF(/i"))  *  *I{x,  t  -f  At) 

(x,  t  +  At)=  LPF(/i"))  *  (x,  t  +  At) 

n  =  n  +  1 

1  (”)  f  r(^~t)  y(T!.— 1)  I  (n)  _  0 

solve  ^  from  -Ip  +^f  '^Ip  —  u 

d(")  =  +  At  vP 

I*'”)  (x,  t  +  At)  -  /(x  -  d^”^ ,  t) 

fin)  ^  2/(-l) 

•  end  iteration  ,  ^ 

v(")  = 

At 


2.6  Overall  conclusion 


Iterative  techniques  were  applied  to  the  estimation  of  motion  in  images.  One  method 
is  based  on  a  direct  formulation  using  the  displaced  frame  difference.  A  method 
which  overcomes  the  scaling  problems  of  the  latter  method  is  based  on  the 
linearization  of  the  displaced  frame  difference,  using  additional  smoothness 
constraints  imposed  on  the  vectorial  velocity.  An  efficiency  improvement  results 
from  the  definition  of  the  normal  image  velocity,  thereby  reducing  the  number  of 
unknowns  with  a  factor  two.  A  continuous  velocity  field  is  an  inherent  property  of 
the  techniques.  Subpixel  displacements  can  be  obtained  with  all  three  methods  and 
multiresolution  schemes  are  also  discussed. 

The  computation  of  gradients,  the  setting  of  different  filter  characteristics  and  the 
optional  use  of  weighting  functions  require  further  study.  Furthermore,  special 
precautions  have  to  be  taken  to  avoid  folding  contributions  due  to  the  periodicity 
associated  with  the  Fourier  or  cosine/sine  transforms. 
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